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Abstract 


Sensitivity  analysis  is  an  integral  part  of  virtually  every  study  of  system  reliability. 


This  paper  describes  a  Monte  Carlo  sampling  plan  for  estimating  this  sensitivity  in  system 
rei :  j  abi^y^t^c^lingel^in  ■component  reliabilities.  The  unique  feature  of  the  approach  is  that 


reliability  tdxhanges  in'eomponent  reliabilities.  The  unique  feature  of  the  approach  is  that 
sample  data  collected  on  K  independent  replications  using  a  specified  component  reliability 
vector  p  are  transformed  by  an  importance  function  into  unbiased  estimates  of  system 
reliability  for  each  component  reliability  vector  q  in  a  set  of  vectors  1 .  Moreover,  this 
importance  function  together  with  available  prior  information  about  the  given  system 
enables  one  to  produce  estimates  that  require  considerably  less  computing  time  to  achieve  a 
specified  accuracy  for  all  1 1  |  reliability  estimates  than  a  set  of  1 1 1  crude  Monte  Carlo 
sampling  experiments  would  require  to  estimate  each  of  the  1 3,  |  system  reliabilities 
separately.  As  the  number  of  components  in  the  system  grows,  the  relative  efficiency 
continues  to  favor  the  proposed  method. 


The  paper  shows  the  intimate  relationship  between  the  proposal  and  the  method  of 
control  variates .  It  next  relates  the  proposal  to  the  estimation  of  coefficient  s  in  a  reliability 
polynomial  and  indicates  how  this  concept  can  be  used  to  improve  computing  efficiency  in 
certain  cases.  It  also  describes  a  procedure  that  determines  the  p  vector,  to  be  used  in  the 
sampling  experiment,  that  minimizes  a  bound  on  the  worst  case  variance.  The  paper  also 
derives  individual  and  Simultaneous  confidence  intervals  that  bold  for  every  fixed  sample 
size  K.  An  example  illustrates  how  the  proposal  works  in  an  s-t  connectedness  problem. 
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Abstract 


Sensitivity  analysis  is  an  integral  part  of  virtually  every  study  of  system  reliability. 
Tots  paper  describes  a  Monte  Carlo  sampling  plan  for  estimating  this  sensitivity  in  system 
reliability  to  changes  in  component  reliabilities.  Tbe  unique  feature  of  the  approach  is  that 
sample  data  collected  on  K  independent  replications  using  a  specified  component  reliability 
vector  p  are  transformed  by  an  importance  function  into  unbiased  estimates  of  system 
reliability  for  each  component  reliability  vector  q  in  a  set  of  vectors  1 .  Moreover,  this 
{importance  'functkM  together  with  available  prior  information  about  the  given  system 
enables  one  to  produce  estimates  that  require  considerably  less  computing  time  to  achieve  a 
specified  accuracy  for  all  1 1 1  reliability  estimates  than  a  set  of  1 1  \  crude  Monte  Carlo 
sampling  experiments  would  require  to  estimate  each  of  the  1 1  \  system  reliabilities 
separately.  As  the  number  of  components  in  the  system  grows,  the  relative  efficiency 
continues  to  favor  the  proposed  method. 

The  paper  shows  the  intimate  relationship  between  tbe  proposal  and  the  method  of 
controi  variates.  It  next  relates  the  proposal  to  tbe  estimation  of  coefficients  in  a  reliability 
polynomial  and  indicates  how  this  concept  can  be  used  to  improve  computing  efficiency  in 
certain  cases.  It  also  describes  a  procedure  that  determines  the  p  vector,  to  be  used  in  the 
sampling  experiment,  that  minimises  a  bound  on  the  worst  case  variance.  The  paper  also 
derives  individual  and  simultaneous  confidence  intervals  that  hold  for  every  fixed  sample 
site  K.  An  example  illustrates  how  the  proposal  works  in  ai  s-t  connectedness  problem. 

Key  Words:  s-t  reliability,  system  reliability,  importance  sampling,  Monte  Carlo 


method,  control  variates 


"")  Sensitivity  analysis,  which  represents  an  integral  part  of  virtually  every  study  of 
system  reliability,  measures  variation  in  this  quantity  in  response  to  changes  in  component 
reliabilities  or  in  system  design.  Replacing  old  components  with  new  ones  with  higher  relia¬ 
bilities  affects  system  reliability.  As  time  elapses,  system  reliability  deteriorates  when  a 
non  replacement  policy  for  component  failures  is  in  force.  Deleting,  adding  or  rearranging 
components  all  affect  system  reliability.  Sampling  variation  in  component  reliability  esti¬ 
mates  induce  sampling  variation  in  the  corresponding  system  reliability  estimate.  Having 
access  to  a  model  that  accurately  predicts  these  changes  in  system  behavior  allows  one  to 
make  considerably  more  well  informed  decisions  for  maintaining  or  enhancing  performance. 


paper  presents  a  method  for  estimating  variation  in  system  reliability  in  res¬ 
ponse  to  variation  in  component  reliabilities.  It  describes  a  Monte  Carlo  sampling  plan  that 
on  each  replication  provides  sample  data  that  contribute  to  the  estimation  of  system  relia¬ 
bility  for  each  of  w  sets  of  distinct  component  reliabilities.  The  sets  may  represent  alter¬ 
native  component  replacement  plans,  deteriorating  component  reliabilities  at  a  succession 
of  time  points  or  extremal  points  of  simultaneous  component  reliability  interval  estimates 
(Fishman  1987).  For  purposes  of  exposition,  we  focus  on  s— t  reliability  but  emphasize  that 
the  concepts  d'acussed  here  also  apply  to  other  definitions  of  system  reliability.  ( I K?ut/sOr* f  * 

To  understand  the  significance  of  this  approach,  we  first  discuss  the  computation  of  7 
s-t  reliability  at  a  point.  Consider  a  system  representable  by  an  undirected  network  C  =*' 
(K&)  where  V denotes  the  set  of  nodes  all  of  which  function  perfectly  and  8  denotes  the  set 
of  edges  each  of  which  fails  randomly  and  independently.  The  concept  of  S-t  reliability 
principally  focuses  on  the  probability  that  at  least  one  path^of  functioning  edges 
(components)  connects  nodes  s  and  t€  7,  and  it  is  this  quantify'that  we  wish  to  compute. 

Since  the  exact  computation  of  9—$.  *di ability  from  a  single  set  of  component 
reliabilities  belongs  to  the  class  of  NF^hard  problems  (Valiant  1979),  attempts  to  shed 
light  oil  this  computation  have  had  to  exploit  special  structure,  rely  on  bounds  or  use  the 
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UoBte  Carlo  method.  The  polynomial  time  Algorithm  of  Agrawai  and  Satyanarayana 
(1964)  far  the  exact  reliability  computation  for  aeriee-parallel  systems  exemplifies  a  highly 
beneficial  me  of  special  structure.  Bounds  such  as  those  of  Eaary  and  Proschan  (1966),  Van 
Syks  cad  Frank  (1972)  and  Ball  and  Provan  (1983)  offer  interval  approximations  to  the 
reliability.  The  Monte  Carlo  method  also  exploits  special  structure  and  bounds.  In 
particular,  the  sampling  plans  in  Van  Slyke  and  Frank  (1972),  Kumamoto,  Tanaka  and 
Inoue  (1977),  Easton  and  Wong  (1980),  Kumamoto,  Tanaka,  Inoue  and  Henley  (1980), 
Karp  and  Luby  (1963)  and  Fishman  (1966)  describe  how  to  derive  statistical  estimates  of 
reliability  that  are  generally  more  accurate  than  a  crude  Monte  Carlo  sampling  experiment 
would  provide  based  on  the  *arre  amount  of  work. 

With  regard  to  the  exact  computation  of  s-t  reliability  for  w>l  sets  of  alternative 
component  reliabilities,  the  corresponding  time  complexity  has  the  same  form  as  that  for  a 
single  reliability  computation  increased  by  the  multiplicative  factor  w.  Also,  whereas  any 
of  the  aforementioned  Monte  Carlo  proposals  allow  one  to  estimate  s-t  reliability  for  each 
of  the  w  points,  regrettably  an  observation  on  a  trial  for  one  point  in  no  way  contributes  to 
estimating  system  reliability  at  the  remaining  w-1  points.  The  present  paper  overcomes 
this  last  inadequacy-  It  describes  a  Monte  Carlo  sampling  plan  that  on  each  trial  generates 
data  that  contribute  to  estimating  all  w  system  reliabilities  simultaneously.  Most 
importantly,  these  estimates,  at  all  w  points,  are  considerably  more  accurate  than 
corresponding  ertimates  that  crude  sampling  can  produce  for  the  nine  amount  of  work. 

The  proposed  method  exploits  importance  sampHuf,  a  technique  that  Kahn  (1950) 
and  Kahn  and  Harris  (1951)  first  described  for  reducing  the  variance  of  a  Monte  Carlo  esti¬ 
mator  nf  a  point.  The  present  account  extends  this  technique  to  reliability  function  estima¬ 
tion  and,  in  particular,  shows  how  one  can  use  knowledge  available  to  the  analyst  before 
experimentation  to  enhance  the  accuracy  of  the  estimated  function  at  all  w  points  for  a 
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ghrcn  sample  sise  K.  Section  1  introduces  relevant  network  nomenclature.  Section  2  then 
describes  system  reliability  estimation  at  a  point  using  crude  Monte  Carlo  sampling,  as  & 
baseline,  ami  a  highly  efficient  alternative  method  described  in  Fishman  (1986).  Section  3 
extends  the  alternative  method  to  simultaneous  estimation  at  all  w  points,  Section  4  shows 
how  to  assess  the  statistical  efficiency  of  the  proposed  method  and  Section  5  shows  the 
intimate  relationship  between  the  proposed  importance  sampling  technique  and  the  method 
•f  control  variates. 

Section  6  offers  an  alternative  interpretation  of  the  estimation  procedure  that 
reveals  its  relationship  to  estimating  coefficients  in  a  reliability  polynomial  and  shows  how 
this  representation  may  save  computation  time  when  a  small  or  moderate  number  of 
components  vary  their  reliabilities.  Section  7  discusses  how  to  perform  the  importance 
sampling  optimaBi  given  the  set  of  component  reliability  vectors  of  interest.  Sections  8 
and  9  derive  individual  and  simultaneous  confidence  intervals.  Section  10  describes 
essential  steps  feu  implementation  and  Section  11  provides  a  comprehensive  example  that 
illustrates  many  of  the  features  of  the  proposed  technique. 

1.  Problem  Setting 

Consider  a  network  G  =  (T,f)  with  node  set  T  and  edge  set  9.  Assume  that 
nodes  function  perfectly  and  that  edges  fail  randomly  and  independently.  Let 
r  *  number  of  distinct  types  of  edges 
q.  *  probability  that  a  node  of  type  i  functions  i=l,...,r 

q  -  (q,,...^) 

y.  «  set  of  edges  of  type  i 

k .  *  |  9i  |  =  number  of  edges  of  type  i 

k  *  (kj,...,^ 

eij  *  jth  edge  of  type  i  i=l,...,r 

x. .  *  1  if  edge  e.j  functions 


i 

X' 


0  otherwise 
k 

‘Si 


E1  x, 


j*i 


cumber  of  functioning  edges  of  type  i 


ik4 


•»x. 


rl’ 


’l*,A 

*1  “  r 

JT«  set  of  all  edge  states  x 
r 


) 


X. 


P(*M)  -  n  ja|[xijqi+(l-x..)(l^qi)]  -^.‘(1^)  1  1 

—  probability  mass  function  (p.m.f.)  of  state  x€  JT 

4{x)  **  1  if  the  system  functions  when  in  state  x 

*  0  otherwise 

g(q).  E  *x)P(x4mi) 

X€Jf 


(1) 


(2) 


■>  probability  that  the  system  functions 

Wo  also  assume  that  G  describes  a  coherent  system.  A  system  of  components  is  coherent  if 
its  structure  function  {#(*)}  is  nondecreasing  and  each  component  is  relevant  (Barlow  and 
Praschaa  1981,  p.  6). 

For  promt  purposes,  the  system  functions  ($(x)**l)  when  at  least  one  minimal  s-t 
path  (s,t€ 9)  exist*  and  fails  (d(x)«0)  when  no  such  path  exists.  Let  1  denote  a  set  of  w 
component  reliability  vectors  of  interest.  Then  the  purpose  of  analysis  is  to  estimate  the 
H  reliability  function  {g(«j),  q€l}. 


2.  Estimation  at  a  Point 

Crude  Monte  Carlo  sampling  offers  a  baseline  against  which  potentially  more 
efficient  sampling  plans  can  be  compared.  Let  denote  K  independent  samples 

drawn  from  (P(x,k,q),  x£J}}.  Then 


(3) 


i«  an  unbiased  estimator  of  g(q)  with 


VM  g|(q)  -  g(<i)(l-g(<|)]/K.  (4) 

To  compute  g^q),  one  performs  K  trials  on  each  of  which  sampling  X  from 
{P(xtk.<l)}  takes  0(|f|)  time  and  determination  of  $X)  takes  0(max(| P|,| *|))  time, 
using  a  depth-first  search  as  described  in  Aho,  Hopcroft  and  Uliman  (1974).  These  are 
worst  case  times.  One  can  also  show  that  the  ’Dean  total  computation  time  has  the  form 


where 


and 


T(%(q))  -  a,  +  KlOj+OjI  »|  +ai(X  M)] 

Oj(  Ji  k,q)  m  E  P(x,k,q)  C(x) 

xeJT 

G(x)  *  expected  search  time  given  the  component  state  vector  x 


The  quantities  aQ  a v  ou  ind  a^( Jt  k,q)  are  machine  dependent. 

All  Monte  Carlo  methods  described  in  the  previously  cited  references  improve  on 
the  variance  (4).  In  i>articular.  the  method  described  in  Kumamoto,  et  al.  (1977)  and 
Fishman  (1966)  Achieves  this  reduction  by  exploiting  bounds  on  the  structure  function 
{#(«)}  it  is  this  approach  that  we  now  describe  and  later  extend  in  Section  3.  We 
follow  the  development  in  Fishman  (1986). 

Suppose  that  there  exist  9-1  binary  functions  {^Jx),  xe  JK)  and  {^  (x),  xeJ}  such 
that 

*L(*)  <  *(*)  <  ^j(x)  V  X-6JT 


Then  the  system  reliability  g(q)  has  lower  and  upper  bounds  gL(q)  and  gu(q),  respectively, 

where 


ie{L,U}. 


Choosing  Bounds 

The  choice  of  bounds  g^(<i)  and  gg(q)  depends  on  the  reliability  computation  under 
consideration.  As  an  example  for  the  s-t  connectedness  problem,  let  ^  denote 


edge-disjoint  minimal  s— t  paths  of  G  and  let  denote  edge— disjoint  minimal  s— t 

cutsets  of  G.  Let 

^M-i-nfi-n  n‘  x  ] 

■  =1L  1  =  1  J  = 1  J 

e. 

ij  i  ■ 


J  r  k  . 

n  i-n  n1  (i-x..) 

=  l[  i  =  i  j  =  l  1J 

e.  .earn* 

ij  i  ■ 


which  are  lower  and  upper  bounds,  respectively,  for  $(x).  Then 


I  r  r  \8.(\9\ 

gL(q)  =  i- n  i-nq.  1  “ 

■  si1-  i  =  l 

and 

j  r  lar.nff  | 

g(J(q)  =  n  i-  n  (l-q  )  1  “ 

■  =  1  [  i  =  l 


are  lower  and  upper  bounds,  respectively,  on  g(q).  One  can  determine  .5®,...,^  in 
0(1 1  y|)  time  using  a  network  flew  algorithm  with  unit  capacities,  as  in  Wagner  (1975, 
p.  954),  and  in  0(|  £|)  time  by  beginning  at  node  s  and  appropriately  labeling 

arcs.  Note  that 

I  <  size  of  the  smallest  minimal  s-t  cutset  in  G 


and 


J  <  size  of  the  smallest  minimal  s-t  path  in  G. 

Moreover,  the  resulting  form  of  {Q(x,k,q)}  in  (5)  enables  one  to  use  Procedure  Q  in 
Fishman  (1986)  to  sample  x  in  0(  |  S\ )  time. 

To  comp  *'?  gK(q)  using  precomputed  bounds  based  on  edge-disjoint  minimal  s-i 
paths  and  cutsets,  one  performs  K  trials  on  each  of  which  sampling  X  from  (Q,x,k,q)} 
occurs  in  0(|£|)  time  using  Procedure  Q  in  Fishman  (1986),  and  determination  of  0(X) 
again  takes  0(max(  |  Y\ ,  |  S\ ))  time.  Also,  mean  total  time  assumes  the  form 


T(g,(q))  «  /J0  +  K|0,+fl,|  y|+a3(^1,k,p)/A(q)] 

where 

=  {xeJT:  *L(x)=0  and  ^(x)=l} 

and  0^...,02  denote  machine  dependent  constants. 

Observe  that 


K(q)  =  K  var  gj(q)/var  gf(q) 


denotes  the  number  of  trials  one  would  have  to  take  with  crude  Monte  Carlo  to  achieve  the 
same  variance  that  arises  in  K  trials  using  {Q(x,k,p)}.  Then  A^q)  = 
T(gf(q)(q))/T(gl(q))  measures  the  efficiency  of  gj(q)  relative  to  g^(q)  and  for  large  K  and 
|  ?|  has  the  approximate  form 


A,(q) * 


;  «,+«,(*  k,q)/lJ|  g(<l)[l~g(q)] 

< -*;M)/A(q)  |  *\  [g^qHtfq)]  [g(q)-gL(q)l 


(9) 


where  Oj(  JT ,k,q)/ 1  £| )  and  o3(J^1,k,q)/ 1  g\  are  bounded  from  above.  A  ratio  greater  than 
unity  favors  the  alternative  sampling  plan.  Experience  (Fishman  1986a)  has  shown  this 
usually  to  be  the  case  by  a  large  margin. 


3.  Estimation  at  a  Set  of  Points 

This  section  extends  the  technique  based  on  bounds  for  a  single  point  to  estimation 
at  a  set  of  w  =  |  £|  points.  In  particular,  it  shows  that  at  least  two  estimators  deserve 
attention  and  later  Section  5  shows  how  a  linear  combination  of  these  estimators  is,  in  fact, 
a  control  variate  estimator.  Let  p=(p.,...,p  )  0<p.<l  for  i=l,...,r  and  let 

a  r  x 


~9- 


R(x£,q,p)  «  P(x^,q)/P(x,k,p) 

=  n  (q./pi)*1  [( i — <ii)/ ( i  Pi )] *i'  <10) 

Lemma  1.  Let  X  be  sampled  from  the  p.m.f.  (Q(x,k,p)}  in  (5).  Then 

EWX)  R(Xk,q,p)]  -  [g(q)~«L(q)]/A(p)  (11) 

and 

E[^(X)  R(X£,q,p)]2  «  c(q,p)(g(q*)-gL(q*)]/ A(p)  (12) 

where 

r  k. 

c(q»p)  —  n  c.1 

i  =  i  1 

c.  =  c(q.,Pi)  *  qJ/pj  +  (1— qi)2/(l— Pi)  **l.-»r 

and 

q*  =  (qj/c^i* -.qJ/CjP,)- 
The  Appendix  contains  the  proof. 

Theorem  1  shows  how  these  properties  relate  to  reliability  estimation. 

Theorem  1.  Let  X  be  sampled  from  the  p.m.f.  (Q(x,k,p)}  in  (5)  and  let 

1&4(x,q,p)  =  gL(q)  +  A(p)  $x)R(x,k,q,p)  (13) 


and 


V\,(^q,P)  =  g„(q)  -  A(p)(l-«(l)]R(x,k,q,p)]  qel. 


(14) 
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Then  for  each  qe  1 

E^(X,q,p)  =  E^(X,q,p)  =  g(q)  (15) 

var  ^/X,q,p)  =  va(q,p)  =  c(q,p)A(p)[g(q*)-gL(q*)]-(g(q)-gL(q)]2  (16) 

var  ^(X,q,p)  =  vb(q,p)  -=  c(q  p) A(p)(gu(q* hlqiH^qHCo)]2  (17) 

and 

cov(V»a(X,q,p),  V>b(X,q,p)]  =  vab(q,p)  =  [g^qHCqJMqHJq)].  (18) 


The  proof  follows  from  Lemma  1. 

Observe  that  the  importance  Junction  R  corrects  the  expectations  (15)  to  the  desired 
value,  thereby  inducing  the  variances  in  (16)  and  (17).  The  implication  is  immediate.  Let 
X*1)  now  denote  the  ith  sample  drawn  from  {Q(x,k,p)}  and  observe  that  one  now  has  two 
potential  estimators  of  g(q),  namely 


with 


-  ,  x  i  1  (0 

gj,(«l,P)  =  k  ,E  ^j<x  >W) 


v" «j,(<l.P)  =  w  VyX,q,p)/K 


(19) 

jf{»,b}.  (20) 


Most  importantly,  note  that  merely  sampling  with  edge  reliabilities  p  enables  one  to 
generate  two  unbiased  estimator  of  the  entire  reliability  function  (g(q),  q€J&}.  By 
contrast,  all  Monte  Carlo  sampling  plans  cited  in  the  introduction  require  |  £\  separate 
experiments  to  estimate  {g(q),  qe.2}. 


4.  Efficiency 

Measuring  the  statistical  efficiency  of  (g^^p),  qei}  and  (gbK(q,p),  q€.2}  as 


**banan  of  {g(q),  q€J*  nib  for  »  more  elaborate  analysis  than  that  for  estimation  at  a 
single  point.  In  particular,  the  sobering  observation  that  c(q,p)  in  (16)  and  (17)  increases 
wpwxBt^iy  with  1 9 1  makes  one  circumspect  about  the  benefit  of  the  proposed  method  as 
the  skee  of  G  grows.  We  now  show  that  this  benefit  is  assured  for  finite  1 1\  and  number  of 
edge  types  r,  provided  that  pel . 

Let  J  «  where  q^  *  (q^,—,^)  end  qjj  is  the  reliability  assigned  to 

anpoasnts  of  type  i  in  the  jth  component  reliability  vector  for  j=l,...,w.  Let  eft  - 
{1,...^}  and 

J9*  *  {i€  p^. j  for  at  least  one  j;  j=l,..., 1 1\ }, 

so  that  |  J9*\  component  reliability  types  vary  in  1 .  Algorithm  A  describes  the  steps  for 

computing  the  estimates  and  provides  the  basis  for  measuring  efficiency.  In  addition  to 

computing  {grf(q,p),  g^p);  qfJ&,  it  computes  {V(ga(q,p)],  V[gM(q,p)];  <£l}  as 

unbiased  estimators  of  {var  g—(q,p),  varg^q.p);  qeJj.  Observe  that  preprocessing  in 

step  1  takes  0(  |  J9*\  |  1) )  time,  postprocessing  in  step  3  takes  0(  |  Jt\ )  time  and,  on  each 

replication,  sampling  in  step  2a  takes  0(|7|)  time,  summation  in  step  2c  takes 

0(  £  k.)  <  0(|  f|)  time,  determination  of  $X)  in  step  2b  takes  0(max|  Y\  ,|  I’D)  and 

i€«^r*  1 

step  2d  takes  0(  |  J9*  ]  |  jf|)  time.  One  can  also  show  that,  the  mean  total  time  for  K 
replications  using  Algorithm  A  has  the  form 

T({ga(q,p),glt(q,p)))  -  «0+u,  |  «** 1 1  +41, 1 1|  +K[oi3+^  I  *|  +a3( J^,,k,p)/A(p) 

+ur4|^r*||g|+a.5  E  kj 

time  where  u>0,...,a>5  denote  machine  dependent  constants  and  /?2  is  identical  with  /?2  in 
T(gf(q)).  To  reduce  numerical  error,  all  computation  in  step  3  should  be  performed  in 
double  precision  arithmetic. 
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Al«sritk-i 

fMPfMt  fa  NtiaU  tk  reliability  f section  {((4),  q«  J}. 

fcprti  Rotvork  f  *  {Kt)i  aarbar  of  type  of  co^oaaats  r;  =  iuber  of  components  of 

type  i  far  i»l,...,r;  saapliaf  distribatioa  {|(*,k,p),  x«J|;  set  of  of*  =  set  of 
ci— t  types  tkafe  nry  ia  &  ;  lover  aad  sppar  boaads  {^(q) ,  ^(q) ;  q«lU{p}} ; 

sad  — bar  of  iadapsad— it  raplicatioaa  I. 

tipis  {fa(q.p)>  IM(q>») .  V[*a(q,p)].  ▼C«bl(q.p)];  “  aabiaaed  estimates  of 

A  A 

(*(q),  s(q).  w  ia(q,p),  w  ^(q.p);  q«$. 

H—fl&ii  • 

■WMa 

1.  Iaitialiaatioa 

••  A(p)  ♦'^(p) 
b.  far  aacb  qti  : 

t(q)-Si(q)*f(q)*fi(q)‘-0. 
for  aacb  itJf*: 

•t(q)  ♦-  la[qi(l-Pi)/P1(l-qi)]  oad  ^(q)  ♦-  la^l-q^/tl-p^] . 

2.  to  aacb  of  I  iadapaadaat  trials: 

a.  Saapla  I  j*l, . . . ,k^  i=t , . . .  ,r  fro-  {Q(x,k,p)}  as  ia  Fiskaaa  (1986) . 

b.  Pataniae  d(I) . 

e.  I  ►  S1 1  . 

|.l  « 

d.  for  q«J> : 

T(q)«-0. 

for  each  ie^f*:  T(q)  «-  T(q)  ♦  k^  (q)  +  X^o^a) . 

T(q)  -  «*p[T(q)] ;  S(q)  -  S(q)  ♦  T(q) ;  V(q)  «-  V(q)  ♦  d(I)  T(q)  ; 

S^q)  -  ^(q)  ♦  T(q)  T(q);  V^q)  -  V^q)  +  P(X)  T(q)  T(q)  . 

3.  Coapatatioa  of  s— ry  statistics 

For  sack  qti : 

*"  *!,(*)  +  V(^)/1- 
««(q>p)  -*u(q)  -  A(p}  pw-*w- 

VC*lI(q.F)]  -  Aa(p)  {[Vt(q)-I[V(,)/I]2}/I(I-l). 

V[«bI(q»f)]  -  ^(pJip^qJ-V^qJl-CCtSCqJ-Viq))/!]2}/!^!). 


tad  of  procadare 


Let  us  now  compare  this  approach  to  estimating  {g(q),  q€  with  the  alternative 

approach  based  on  the  |  Jj  point  estimates  {g^  \(q),  <je  J}  using  (3),  where  one  chooses 

the  sample  siaes  (K(q,p),  q€l}  to  achieve  equal  vanances.  That  is, 


«KV)(q)  *  S«?]/K(q,p) 


where 


Observe  that 


K(q,p)  *  K  A(q,p) 


w  v  *(q)(l-g(q)l 

*  min - var  0.(< 


A(p,p)  *  g(p)[l-g(p)]/[gc(?)-g(p)][g(p)-gL(p)] 


and,  except  in  special  cases,  for  any  odgs  type  i€eV* 


lim  A(q,p)*0 
k^o# 


for  q^p. 


This  last  limit  follows  from  the  growth  of  c(q,p)  with  k.. 


A(p)  *  £  A(q,p). 
q6JE 

and  observe  tuai, 


lim  A(p)  *  A(p,p). 
k.  -»oo 


Observe  that  the  time  ratio 
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\(J>P) 


T({ga(q,p),gbI(q,p)}) 


(24) 


where 


measure*  the  effidency  of  the  proposed  method  relative  to  using  crude  Monte  Carlo 

templing  with  (3)  |  Jj  time*  to  obtain  estimate*  with  equal  variances  var  ^(q)  = 

min  var  g »(q,p)  for  qe  1 .  As  k.  increases,  (24)  assumes  the  form 
j  *  1 

E  J  «2+aj(  jr»k><i)  /ki)A(q  ,p) 

A  (  J  ,p)  «  - .  (25) 

(  *^1 » k>p)/A(p)ki+w5 

Practice  indicates  that  u^«0r  Then  provided  pc  1 ,  (25)  generally  satisfies  ,p)  > 
Aj(p)  for  large  k.  (ic*r*),  implying  that  Algorithm  A  is  at  least  as  efficient  as  (6)  at  a 
single  point.  As  the  example  in  Section  11  shows,  the  realized  efficiency  can  be 
considerably  greater. 

5.  Tim  Optimal  Estimator 

The  representations  of  grf(q,p)  and  gM(q,p)  in  (19)  suggest  that  these  quantities 
conceptually  are  alternative  forms  of  a  more  comprehensive  estimator.  Theorem  2  confirms 
this  observation. 

Theorem  2.  Let  X  denote  a  sample  from  (Q(x,k,p)}  and  define 

tK*,q,p,©)  =  e^(x,q,p)  +  (i-0)  ^(x^p)  -  »<e<®.  (26) 
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Tkea  for  fixed  p* 

i.  Ed(X,q,p,0)  -  g(q) 

»•  ©*(«)  -  (vfc(<fcPH4b(<ltP)]/[v%(q,p)+vb(q,p)-2v4b(q,p)) 

-  i/{i+l\(«.pWBb(q»p)!/(vb(q,pHab(q,p)l} 

minimises  var  ^X,q,p,8) 

Ui.  var  tfX,q,p,0*(q,p))  «  (va(q,p)vb(q,p)-vJb(q,p)]/[c(q,p)A(p)A(q*)-A2(q)]. 

The  proof  follows  directly  from  the  minimisation  of  var  #X,q,p,6)  with  respect  to  8. 
Therefore,  among  estimators  of  the  fonn 

i,(q,p,o)  *  ®srf(«i,p)  +  (i-9)  ^(q.p),  (27) 

(q,p))  has  minimal  variance. 

Observe  that 

©*(q»p)  *  1  if  va(q>p)  =  v^fop)  t  vb(q,p) 

=  0  if  vb(q,p)  =  vab(q,p)  i  vjq,p), 

revealing  that  g^(q,p)  is  optimal  it 

vsr  ga(q,p)  =  var  gf(q) 

and  gtf(q,p)  is  optimal  if 


v«  gM(q,p)  =  var  gr(q) 
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whw  g^q)  and  var  |g(f)  are  defined  in  (6)  and  (?)  respectively. 

Writing  (26)  in  the  alternative  form 

#x,q,p,0)  *  ^(x,q,p)  +  6[R(x,k,q>p))-A(q)]  (28) 

tmta  that  . R<X(I)M,P)  act  as  control  m notes  with  known  mean  A(q) 

and,  by  appropriate  robetitntion,  that 

**  d(X,q,p,0*(q,p))  *  var  ^)(X,q,p){l-corr2(^>(X,q,p),R(X^,q,p)]}, 

where  corr(A3)  denotes  the  coefficient  of  correlation  between  A  and  B.  Note  that  this 
variance  diminishes  as  the  correlation  between  ^(X,q,p)  and  R(X,k,q,p)  increases  in 

magnitude. 

In  practice,  {0*(q,p),  q E  J}  is  unknown  but  can  be  estimated  unbiaaedly  for  q^p  by 
6*(«i,P)  -  X{V[gbt(q,p)}-ZI(q,p)}/[c(q,p)A(p)A(qVA\q)]  (29) 

where 

V^p)  *  i%(o>-8bI(<i>p)]i8rf(q.pHl(<i)]/(K-i) 

V(gM^q,p)]  is  computed  in  step  3  of  Algorithm  A.  Although  one  may  incline  to  use  the 
estimator 


«f(>q»p^  \n.ir))  ~  fea(q,p)  +  [i-Q*(q,p)]gh^q,p) 


(30) 


for  g(q),  one  ouickly  sees  that 
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®ig(*P.6*(%p))  -  g(q)  +  cov  (S*(q,p),  g^p)) 

-  cov  (d*(q,p),  g^q.p)]- 

Therefore,  jf(q^,4*(q,p))  generally  it  biased.  While  this  bias  diminishes  as  K  increases, 
the  rate  of  dissipation  depends  on  the  particular  system  under  study.  As  a  result,  no 
general  statement  it  possible  regarding  how  large  K  must  be  in  order  to  treat  bias  as 
lackfontal.  Because  of  this  limitation,  the  remainder  of  this  paper  focuses  on  choosing 
between  grf(q,p)  and  g^ (q,p).  We  return  to  the  issue  of  choosing  a  point  estimator  for 
g(q)  in  the  example  of  Section  11. 


6.  An  Alternative  Representation 

This  section  describes  an  alternative  representation  for  g(q)  that  offers 
considerable  conceptual  value  for  function  estimation.  Let 


and 


K  . 

JT(s  )  *  {*eJT|  Elx  *  g  i~l,...,r} 

1  r  j=l  1 

=  E  d(x) 

x€jr(slt...,ir) 

3=  the  number  of  possible  ways  that  the  system  can 
function  (d(x)=l)  when  *lt...,zr  components  of  types 
l,...,r,  respectively,  function  and 
components  of  types  l,...,r,  respectively,  fail. 


(31) 


(32) 


Then  one  can  write 


*(«)  =  OVV-  MrM*, . \> 

*,=01  *  =0  r 


(33) 


0,1 . k;  0<0<1,  k=l,2,... 


f.(M)  -  (J)  ^(l-q)^  «- 

and 

tt(*1.»M*r)  -  M(*lt...4r)/  n  £). 

i*l  i 

Obwiw  that  since  we  are  working  with  a  coherent  system 
u(0,...,0)  ■  0 

u(kj,...tkr)  *  1 , 

so  that  0<s. <kit  1*1, ...,r}  is  t  multivariate  distribution  function. 

If  one  were  to  perform  crude  Monte  Carlo  sampling,  then  the  estimation  of  g(q) 
would  be  equivalent  to  eetimating  the  coefficients  {u^,...,*^}.  In  the  case  of  r=l,  M(Zj) 
denotee  the  number  of  connecting  cutsets  of  G  when  z^  arcs  function  and  arcs  fail. 
Ftehman  (1987)  describes  a  method  of  estimating  M(z1)  in  this  special  case.  As  we  now 
show,  the  present  proposal  corresponds  to  the  implicit  estimation  of  analogous  quantities. 
Let 


k. 

n  (tl) 

_  4  m  - 


i  =  1 


1  E  4  (x)  je{L,U) 

*€**(*, . zr )  J 


(34) 


*  number  of  replications  on  which  the  system  (35) 

functions  when  zt,...,*r  components  of  types  l,...,r, 
respectively,  function  and  k1-z1,...,kr-«r  components 
of  types  l,...,r,  respectively,  fail 


and 


(36) 


-  number  of  replications  on  which  the  system 

fails  when  ij,...,s  components  of  types  l,...,r, 
respectively,  function  and  -*  components 

of  types  l,...,r,  respectively,  fail. 


Then  grf(q»p)  uid  gu(q,p)  in  (6)  have  the  equivalent  forms 


g„(q.P)  -  +  A(p)  £*...  ErR''(»,k,q,p)  . . . 

*j-0  *r“0 


(37) 


and 


kt  k  * 

8M(<l.P)  *  g0(q)  -  A(p)  E1...  ErR  (z,k,q,p)  K  (z  .,z  )/K 

*,=0  *  =0  ft  l  r 

1  r 


(38) 


whan 


R*Ww)  -  n  (qj/p.r  |(l-qi)/(l-pi)]kr*i 

i.l  *  1 

EK4(Sj,...,*r)/K  =  u^Sj,...^^)  =  u(Zj,...,*r)  —  u^(*j,...,zr) 


and 


~  Ub(Zl’"’’Zr)  u(Zj,...,Zr). 

Therefore,  using  g^fop)  is  equivalent  to  estimating  the  coefficients  {ua(zlt...,zr)} 
implicitly  and  using  g^qp)  is  equivalent  to  estimating  the  coefficients  {ub(zr...,zr)} 
implicitly. 

Expressions  (19),  on  the  one  hand,  and  (37)  and  (38),  on  the  other,  have  beneficial 
and  limiting  features.  If  one  uses  Algorithm  A,  then  the  sample  reliability  functions 
8bj(q»P)i  available  for  study  have  ordinates  only  at  the  points  in  2, 
specified  in  the  sampling  experiment.  However,  if  one  alternatively  records  (Ka(zr...,zr), 
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then  (37)  and  (38)  enable  one  to  construct  the  sample  functions  for  any  q  in 
[0,l]r  which  may  be  of  interest  at  any  time  after  the  sampling  experiment  terminates.  Also, 
aa  Section  9  shows  shortly,  for  a  given  network  G,  the  widths  of  confidence  intervals  for 
(g(q),  J)  based  on  either  {ga(q),  q€^  or  (g^q),  q€  J)  as  in  (19)  increase  with  |  £\ 
whereas  confidence  intervals  based  on  {K  (zir..,zr),  Kb(zlt...,zr)}  have  widths  independent 
of|J|. 

This  alternative  approach  requires  a  space  of  II  k.  counters  to  accumulate 

i€*ST*  1 

(Ka(*)}  end  a  like  space  to  accumulate  {Kb(- )}.  Moreover,  a  modified  Algorithm  A  based 

on  this  approach  would  replace  the  time  components  0(K  |  <4f*|  1 1\ )  and  0(K  E  k.)  in 

i  €<ar*  1 

T({gll(q,p),  g,*(q,P)})  by  0(|  J|  n  k.)  <  0(|  *|/|  |)'  thus  eliminating  their 

dependence  on  the  sample  size  K.  If  II  k.  is  large,  this  may  limit  the  extent  to  which 

ie<ar* 1 

one  can  store  the  sample  sums  (Kt(-),  Kfe(  • )} .  However,  when  they  are  storable,  their 
availability  offers  considerable  post-experimental  discretion  for  computing  quantities  of 
interest. 

7.  Choosing  the  Sampling  Vector  p 

The  forms  of  the  variances  (16)  and  (17)  cle^ly  indicate  that  the  choice  of  p  affects 

the  statistical  accuracies  of  g^^p)  and  gM(q,p).  While  no  unequivocal  rule  exists  for 

choosing  p,  minimizing  max  [min  var  g._(q,p)]  is  one  reasonable  objective.  Unfortu- 

q€J  j€{*,b} 

nately,  the  unknown  variances  render  this  minimization  impossible.  An  immediate  alter¬ 
native  uses  the  upper  bound 

c(q,p)  A(p)  A(q')  >  max  var  giK(q,p) 

j  «{a.b}  JK 

and  finds,  by  grid  search,  the  p  that  minimizes  A(p)  max  c(q,p)  A(q*).  However,  a 

q6  J 
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considerably  better  bound  h(q,p)  exists  when  either  p<q  or  p>q  so  that  sampling  with  the  p 

that  minimizes  max  h(q,p)  can  produce  considerably  better  worst  case  results. 
qe.2 

This  section  derives  h(q,p)  explicitly  for  var  ^b(X,q,p).  It  makes  use  of  the 
observation  that  for  a  coherent  system  p<q  for  all  qe  1  implies  g(p)  <  g(q)  and  gj(p)  <  gj(q) 
for  je{L,U}  and  p>q  for  all  qe£  impUes  g(p)  >  g(q)  and  g^p)  >  g^q)  for  j{L,U}.  A 
completely  analogous  approach  holds  for  var  g^(qp).  Also,  the  Appendix  extends  the 
analysis  (Theorem  4  and  5)  to  cases  in  which  the  coefficients  of  variation 

71;j(q>P)  =  (var  0j(X,q,p)]Vg(q) 
and 

72j(Q»P)  =  (var  ^.j(X,q,p)]V(l-g(q)]  je{a,b} 

are  the  criteria  of  accuracy. 

Lemma  2.  Define 


hJ(z)  =  h^Zjq  p)  =  c(q,p)A(p)A(q*)  -  [g^qM2 
and  (^9) 

h2(z)  =  h2(z,q,p)  =  c(q,p)A(p)[gu(q*)— z]  -  [gw(q)-z]2 

—  00<z<00. 


Then  for  either  p<q  or  p>q 


-22- 


var  A(X,q,p)  <  h(q,p)  *  max'  max  *  h  (z),  *  max  h  (z)  (40a) 

^  [slW  <  z<gL(q )  1  *,,(«  )<*<g„(<i) 2 

»f  gJqKgJq’^g^q) 

as  max  h,(z)  otherwise.  (40b) 

*,.(<!)<  *<g^q)  1 

The  Appendix  contains  the  proof. 

Theorem  3.  For  ht  and  h2  as  defined  in  (39),  p>q  or  p>q,  and 

2*  =  go(q)  “  c(q,p)A(p)/2, 

var  ^(X^P)  <  Mq,P)  =  maxPl1(%,(q*))»  ymax(z*,gL(q*)))] 

if  «l(q)lgl(<r)<gp(q)  (41) 

=  h^g^q))  otherwise. 

See  the  Appendix  for  the  proof.  To  derive  the  analogous  upper  bound  for  var  y>a(X,q,p), 

*  *  * 

one  replaces  z  by  1-z,  gL(q)  by  l-g„(q),  gD(q)  by  l-gL(q),  gL(q  )  by  l-g„(q  )  and  gy(q  ) 

by  l-gL(q*)  everywhere  in  (38),  (39)  and  (40). 

Recall  from  Section  4  that  choosing  p  from  1  is  beneficial  from  the  viewpoint  of 

efficiency  as  any  of  the  k.  grows.  Then  one  can  compute  max  h(q,p)  by  enumeration 

1  qe.2 

for  every  p  in  £  and  select  the  p  that  minimizes  the  maxima.  In  total  |  Jl  |  points  are 
evaluated.  As  the  example  in  Section  11  shows,  this  method  of  choosing  p  can  lead  to 
significant  improvements  in  statistical  efficiency. 

8.  Individual  Confidence  Intervals 

Although  the  distribution  of  (gjj(q,p)-g(q)]/[var  g^(q,p)]'^  converges  to  the 


standard  normal  distribution  as  K  -»  <»,  this  result,  at  best,  can  only  lead  to  a  rough 
confidence  interval  for  g(q).  To  avoid  the  errors  of  approximation  inherent  in  the  normal 
approach  to  confidence  intervals,  we  use  an  alternative  technique. 

Theorem  6.  Let 


Jf  =  {x6JT:  ^L(x)=0,  $x)=l}, 

RJ (Kq.p)  =  max  R(x,k,q,p), 

*  x€JT 


Yj(q,p)  =  feaI(q,p)-gL(q)]/A(p)Ra(k,q,p), 
m(z,w)  =  z  log(w/z)  +  (1-z)  log  [( 1 — a;) / ( 1 — z)]  0<z,uKl, 


let  <i^z,$/2,K)  denote  the  solution  to  m(z,«)  —  ^  ln(£/2)  for  fixed  ze(0,l]  and  6e(0,l),  and  let 


u>  (z,6/2,K)  =  w(z,$/2,K) 


=  0 


if  0<z<l 


otherwise. 


(42) 


Then,  the  interval 


(gL(q)+A(p)Ra(k,q,p)w(YI(q,P)^/2,K),gL(q)+A(p)Ra(k,q,p)w(l-YK(q,p)^/2,K))  (43) 


covers  g(q)  with  probability  >  1  —  6. 


Theorem  7.  Let 


Jth  =  {xeJT:  #x)= 0,  ^(x)=l}, 

RJk,q,p)  =  max  R(x,k,q,p), 
xejrb 
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and 

Then  the  interval 

(gu(q)-A(p)Rb(k>q»P)w*  (l-ZI(q,P)»^/2,K),gu(q)-A(p)Rb(k,q,p)u;’>  (ZI(q,p),5/2,K))  (44) 

covers  g(q)  with  probability  >1-6. 

Proof  of  Theorems  6  and  7.  Inspection  of  (13)  and  (14)  makes  clear  that 
pr[gL(q)  <  ^a(X,q,p)  <  gL(q)+  A(p)Ra(k,q,p)]  =  1 
and  (45) 

Prlgu(q)-A(p)Rb(k,q,p)  <  ^(X,q,p)  <  g0(q)|  =  1. 

The  resulting  confidence  intervals  follow  from  Theorem  1  in  Fishman  (1988)» 

Although  these  intervals  generally  are  wider  than  the  corresponding  normal 
confidence  intervals  would  be  for  given  K  and  6,  they  are  fine  of  the  error  of  approximation 
inherent  in  normal  intervals. 

To  use  these  intervals  in  practice,  one  needs  to  know  {R^(k,q,p),Rb(k,q,  p);  qM}- 
Theorems  8  and  9  formulate  mathematical  programs  aimed  at  computing  these  quantities. 
Since  experience  with  several  networks  for  the  s-t  connectedness  problem  with  £  = 
{qj<  ...  <q^}  has  shown  that  p  -  qx  usually  minimizes  the  worst  case  bound  (42),  we  focus 
on  the  case  p<q, 

Theorem  8.  Let  35  denote  the  set  of  all  minimal  s-t  cutsets  of  smallest  cardinality,  let 


=  log  lqi(l-pi)/pi(l-<li)] 

ie^r* 

and  assume  q.  >  p.  for  V  i eJt*.  Then 

n  ri 

k.-E  z. 

1  jeJ  1 

R  (M*p)  =  n  (qi/pi)  3 

*  i  €<ar*  1  1 

E  z* 
j e8.  3 

t(Mi)/(l-Pi)]  1 

(46) 

where  *  solves  the  integer  program 

min  E  a.  E  z. 

*  i€<**  1  j€*.  3 

(47a) 

subject  to 

E  z.  >  1 
j€*  J  ' 

V  9ZJH 

a 

(47b) 

and 

1 

*> 

VI 

* 

W  w 

•n 

v 

a 

(47c) 

v  je  *(,**). 

(47d) 

The  Appendix  contains  the  proof. 


Theorem  9.  Let  ^denote  the  set  of  all  minimal  s-t  paths  of  smallest  cardinality,  let 


i=l,...,J:  *  C  *(**)} 

■\  =  {ss?-.  j>c  »(<»*),  |^|  =  i^i) 


qA  >  p4  for  V  i€**.  Then 


_  *  _  * 

k.— E  z.  E  z. 

1  3  jear  3 

R^M*)  =.n  (Vpi)  [(i^Qi)/!1— Pi)l 

X  C 


where  s  solves  the  integer  program 


min  E  a.  E  z. 
s  i€<**  1  jef.  1 


(49a) 


subject  to 


E  z.  <  |*|  -1 
j€*  3 

(49b) 

E  z.  >  1 
je*  3 

VJ»e^b 

(49c) 

Z.€{0,1} 


Vje  *(<**).  (49d) 


The  proof  follows  analogously  to  that  for  Theorem  8. 

Recall  that  since  are  edge  disjoint,  I  <  I  =  the  size  of  the  minimal  s— t 

cutset  of  smallest  cardinality.  Therefore,  if  |  Jt  |  <  I*,  then  |  JT  |  =  0  so  that  the 
constraints  in  (47c)  vanish  and 


RJ (M»p)  =  [  n  (q./p.)]  exp(  -  E  min  a.). 
*  i€<**  1  1  9lA t  i$9  1 


(50) 


The  case  of  |  JCj  =  I  requires  more  detail.  If  the  minimal  s-t  cutsets  in  are 

edge-disjoint,  then  (47)  has  the  form  of  a  transportation  problem  with  <  J*  =  the  size 

of  the  minimal  s-t  path  of  smallest  cardinality  and  can  be  solved  using  a  special  purpose 

algorithm  as  in  Dantzig  (1963,  p.  308).  If  the  cutsets  are  not  edge-disjoint,  JT  potentially 

can  have  an  exponential  number  of  members,  limiting  one's  capacity  to  enumerate  them 

all.  This  possibility  suggests  an  iterative  approach. 

Suppose  one  begins  by  relaxing  (47c).  This  gives  the  candidate  solution  (50).  If  the 

set  of  arcs  chosen  there  do  not  form  a  minimal  s-t  cutset  in  3,  then  the  problem  is  solved. 

If  they  do  form  a  cutset  **,  then  one  activates  the  corresponding  constraint.  Let  i*  denote 

the  edge  in  tf*  with  the  largest  a..  Then  [  H  (q./p.)]  exp(  —  E  £  min  a.) 

1  1  1  ie^»\{i*}  1 

solves  the  problem  provided  that  the  selected  edges  do  not  form  a  minimal  s-t  cutset  in  3 

.  If  they  do  form  a  cutset,  then  continued  iteration  becomes  more  complicated  and  one 

may  elect  to  drop  one  of  the  edge-disjoint  paths  in  from  the  lower  bound  gL(q) 

thereby  reducing  the  size  of  A  and  making 

R  (M,p)  =  l  n  (q./p.)]  exp(  -  E  £  min  a.) 

*  ie^r*  1  1  9zA\9niz9  1 

a'  o 

the  solution. 

The  solution  to  (49)  proceeds  in  an  analogous  manner.  If  Ab  <  J*,  then 

Rv(Mj>)  =  [  n  (q./p.)]  exp(  -  £  £  min  a  ).  (51) 

b  ie^*  1  1  1 

If  |  A^  |  =  J  ,  then  one  can  either  drop  a  cutset  in  A^  from  the  upper  bound  gu(q)  and 
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R^(k»<u>)  =  [  n  (qj/Pi)]  exp(  - 

X  t  ^ 


£  £ 
*E*^\*o  ie* 


min  a.) 

i' 


as  the  solution,  or  again  proceed  iteratively.  With  (49c)  relaxed,  (51)  is  the  candidate 

solution.  If  the  set  of  selected  edges  do  not  form  a  path  in  then  (51)  is  the  minimum. 

If  they  do  form  a  path  9  *  with  edge  i*  giving  the  largest  a^  then  the  solution 

[  II  (q./p.)]  exp(  -  £  £  min  a.)  needs  to  be  checked,  etc.  One  anticipates 

i€^r*  1  1  i€tr\{i*}  1 

that  choosing  edge-disjoint  paths  and  cutsets  such  that  ^  and  ^ 

are  empty  generally  will  have  small  effect  on  the  bounds  gL(q)  and  gu(q)  for  large 
networks. 

9.  Simultaneous  Confidence  Intervals 

Although  each  confidence  interval  in  Section  8  holds  with  probability  >  1-5,  the 
Jint  confidence  intervals  for  (g(q),  q€  hold  simultaneously  only  with  probability  > 
1-|  J|5.  This  result  follows  from  a  Bonferroni  inequality.  See  Miller  (1981,  p.  8).  To 
restore  the  joint  confidence  level  to  1-5,  one  replaces  log  (5/2)  by  log  (5/2 1  £\ }  in  (43)  and 
(44)  and  determines  the  corresponding  solutions.  The  effect  of  this  substitution  is  to 
increase  the  constant  of  proportionality  in  the  approximate  interval  widths  from 
[•  *(2/«)]*  to  [2log(2|J(/«)l*  (see  Fishman  1986).  For  5=.01  and  |  4)  =20  one  has 
.  3(2|J|/5)/log(2/5)]i  =  1.25.  For  5=.01  and  |J|=100,  it  is  1.37  and  for  5=  .01  and 
|  Jf|=1000  it  is  1.52.  Moreover,  if  1  denotes  a  continuous  region  in  the  |  if  |— dimensional 
hypercube  (0,1)  I  ^ ,  then  the  resulting  confidence  intervals  have  infinite  widths  and  are 
therefore  useless. 

An  alternative  approach  derives  simultaneous  confidence  intervals  for  (g(q),  q€.£} 
using  the  representation  of  g(q)  in  (44).  In  particular,  it  implicitly  finds  simultaneous 


confidence  intervals  lor  the  coefficients  {u. of  which  there  are  N  <  n  k.  in 

3  1  ^  "  i€<*r*  1 

(94).  For  convenience  of  notation  we  take  |  Jt  *|  =  r  but  note  the  relatively 
straightforward  adjustment  for  |  *T*!  <  r.  Let  *  =  (*j  and  recall  the  definitions  of 
Kfc(«)  and  Kb(«)  in  (35)  and  (36).  Then  {(w(K.(*)/K,4/2N,K),  w*  (1-K.(*)/K,4/2N,K); 
V «},  where  «*( •,*,•)  is  defined  in  (42),  provide  confidence  intervals  for  (u^i)}  that  hold 
simultaneously  with  probability  >  1-4. 

Observe  that  all  coefficients  u^(z)  are  nonnegative  and  that  {u>* (K^.(s)/K,4/2N,K), 
u  (1— Kj(s)/K,4/2N,K}  are  independent  of  q.  Therefore,  for  all  q€4 

(g(q)-A(p)  E1...  EtR*(z,k,qJp)w(K  (*)/M/2N,K), 

V°  V° 

(52) 

g.(q)+  A(p)  5?...  EfR*(«A.<l.p)w*(l-K  (i)/K,i/2N,K)) 

s.=0  *  =0 
1  r 

simultaneously  covers  (g(q),  qeij  with  probability  >  1-4  and  likewise 


k,  k  *  * 

(g^l)+A(p)  r...  E'R  (*,k,q,p)w  (l-Kl(«)/K,i/2N,K), 


i, =0  m  =0 
1  r 


g_(q)-A(p)  E1...  ErR'(.,k,q,p)u’(K.(«)/K,«/2N,K)) 


*r°  V° 


simultaneously  covers  (g(q),  qGi}  with  probability  >  1-4. 

The  most  desirable  feature  of  this  alternative  approach  is  that  the  resulting 
intervals  are  unaffected  in  width  or  confidence  level  by  the  size  of  1 .  However,  since  the 
number  of  quantities  K^(x)  to  be  collected  is  0(  II ^k.),  this  alternative  approach 

becomes  less  feasible  to  implement  as  the  k.  and  r  increase. 


-30- 


10.  Stn  for  InDimottatioo 

To  implement  the  proposed  sampling  plan  to  estimate  reliability  for  s— t 
connectedness,  one  proceeds  as  follows: 

1.  Determine  a  set  of  edge-disjoint  minimal  s— t  paths 

2.  Determine  a  set  of  edge-disjoint  minimal  8— t  cutsets  ffj,...,  tfr 

3.  Compete  (g^q),  g^q);  q6i}. 

4.  Determine  a  sampling  vector  p  from  las  in  Section  7. 

5.  Using  Algorithm  A,  perform  K  independent  replications. 

6.  For  each  qe  A:  compute  RJk^p)  if  V^fop)]  >  V[gbI(q,p)];  otherwise  compute 
RJk,q,p)  (Section  8). 

7.  Using  the  bounds  (RJk,q,p);  qei}  or  {Rb(k,q,p);  qeJj  in  step  6,  compute 
individual  or  simultaneous  confidence  intervals  for  (g(q),  qei}  (Sections  8  and  9). 

Although  these  steps  require  more  work  than  crude  Monte  Carlo  function  estimation 
does,  one  can  develop  computer  programs  with  sufficient  generality  to  compute  all 
quantities  in  steps  1  through  7  for  many  different  network  designs  Reusing  the  programs 
enables  one  to  distribute  the  fixed  cost  of  their  development  over  all  such  network,  making 
the  cost  per  network  incidental. 


11.  Rsampfo 

An  analysis  of  the  network  in  Fig.  1  illustrates  the  proposed  method.  The  network 
has  30  edges  and  20  nodes.  Also,  the  example  assumes  r=l  so  that  all  edges  have  identical 
reliabilities,  allowing  us  to  write  q=q.  Note  that  any  other  specification  with  r>l  can  also 


Insert  Fig.  1  about  here. 


be  accommodated  easily.  The  objective  is  to  estimate  (g(q),  q=.80+.01(i-l)  i=l,...,20} 
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whgi  g(q)  m  probability  that  nodes  s=l  and  t=20  are  connected  when  edge  reliabilities  are 
q.  for  sampling,  we  use  p=p,  again  merely  as  a  convenience.  The  selected  edge— disjoint 


paths  and  cutsets  are 

JJ  -  {3,9,18,27,28} 
»  {2,7,15,24,30} 
-  {28,29,30} 

*4  «  {4, 5, 6, 7, 8, 9} 


i*  =  {1,5,12,21,29} 

=  {1,2,3} 

*3  =  {11,12,14,15,17,18} 
*5  =  {19,21,22,24,25,27}. 


As  a  preliminary  step,  Table  1  shows  the  worst  case  upper  bound  on  var  ^b(X,q,p) 
as  given  in  (41).  Observe  that  the  choice  p=.80  minimises  this  worst  case  bound  and  it  is 

Insert  Tables  1,2,  and  3  about  here. 

this  component  reliability  that  we  use  for  sampling.  A  parallel  analysis  for  var  ^(X,q,p) 
also  chose  p=.80. 

Table  2  compares  the  estimates  of  var  g^Cq.p)  and  var  gy^p)  for  a  sample  size  K 

*  * 

=  1048576  and  shows  the  estimated  control  variate  coefficient  6  (q,p),  as  in  (29).  These 
results  strongly  favor  relying  on  gbt(q,p),  if  the  choice  is  between  this  quantity  and 
i  (q,p).  Table  3  shows  the  resulting  estimates  in  col.  1  along  with  variance  estimates  in 
col.  5  and  individual  99%  confidence  intervals  in  cols.  6-8.  In  contrast  to  the  exact  results 
in  col.  3  which  took  slightly  more  than  one  hour  each  to  compute,  all  results  in  cols.  1,2,4 
and  5  took  72.7  minutes  in  total,  or  4.16  milliseconds  per  replication.  Computation  of  the 
confidence  intervals  took  incidental  time.  Whereas  the  calculated  exact  results  in  col.  3 
were  accurate  to  sixteen  significant  digits  (reduced  to  four  digits  here  for  comparative 
purposes),  the  confidence  intervals  suggest  an  accuracy  to  two  significant  digits  at  the  .99 
level.  If  two  significant  digits  is  acceptable  for  purposes  of  analysis,  then  the  Monte  Carlo 
approach  clearly  prevails. 
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Table  4  shows  the  effect  of  sampling  at  an  arbitrary  point  p=.90  rather  than  at 
p».80.  Although  sampling  at  p=.90  does  produce  better  results  for  q>p,  the  deficiency  of 
sampling  at  p*.80  in  this  interval  is  considerably  less  than  the  corresponding  deficiency  for 
sampling  with  p«.90  for  .80<q<.89. 


Insert  Fig.  2  about  here. 

Figure  2  displays  several  variance  ratios  that  reveal  how  {^^,.80)}  performs 
compared  to  the  crude  estimator  {g^q)}  in  (3),  the  estimator  {^(q)}  in  (6),  and  the 
approximately  optimal  estimator  {^(^,©*^,.80))}  in  (30).  First,  note  that  {gM(q,p)} 
performs  almost  as  well  as  {g|(qtp,©  (q,.80))}.  Second,  observe  that  uniformly  superior 
ratio  for  (g^q.p)}  when  compared  to  (g^q)}.  In  particular,  note  that  these  ratios  exceed 
100  for  q>.95. 

We  now  turn  to  the  efficiency  measure  (22).  Since  V[ga(q,.80)]  >  V[gM(q,.80)]  for 
all  q€4  Fig.  2  makes  clear  that  A^p)  >  105,  indicating  the  clear  superiority  of 
{^(Qv^O)}  over  the  crude  estimator  (g^q)}  in  (3). 
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.76 

.83 

.52 

.89 

.1561D+03 

.78 

.84 

.54 

.89 

.6S67D+02 

.80 

.84 

.56 

.89 

.2953D+02 

.82 

.80 

.58 

.88 

.1363D+02 

.84 

.80 

.60 

.88 

.6519D+01 

.86 

.80 

.62 

.88 

.3222D+01 

.88 

.80 

.64 

.88 

.1545D  f  01 

.90 

.80 

.66 

.87 

.8754D+00 

.92 

.80 

.68 

.87 

.4824D+00 

.92 

.80 

.70 

.86 

.2790D+00 

.96 

.80 

.72 

.84 

.1709D+00 

.98 

.80 

.74 

.83 

.1116D+00 

^q.  =  q  in  2  at  which  h(q,p)  achieves  its  maximum  for 


.4985D-01 

.3234D-01 

.8818D-01 

.1090D+00 

.1642D+00 

.3255D+00 

.6961D+00 

.5410D+01 

.9111D+02 

.1400D+05 

.7866D+09 

specified  p. 
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Table  2 

Comparison  of  V[gaI(q,p)]  and  V[gbK(q,p)] 
(p=.80,  K= 1048576) 


q 

V(ga(q,p)] 

^(Sy[(q*p)] 

©*(q,p) 

.80 

1.00 

— 

.81 

1.58 

-.4941 

.82 

3.04 

-.1799 

.83 

6.07 

-.0852 

.84 

11.61 

-.0440 

.85 

20.97 

-.0233 

.86 

35.92 

-.0122 

.87 

58.78 

-.0060 

.88 

92.67 

-.0026 

.89 

142.10 

-.0008 

q 

V^^p)] 

Mibifa’P)] 

©*(q>p) 

.90 

214.0 

.0001 

.91 

319.0 

.0005 

.92 

475.5 

.0006 

.93 

713.8 

.0006 

.94 

1097.0 

.0005 

.95 

1755.0 

.0003 

.96 

2995.0 

.0002 

.97 

5727.0 

.0001 

.98 

13665.0 

.0001 

.99 

57108.0 

.0000 

Reliability  Estimation 


Provided  by  J.S.  Provan  using  an  algorithm  based  on  cutset  enumeration. 
Computed  as  in  Algorithm  A 


Fig.  1  Network 


All  component  reliabilities  are  identical. 

p  =  .80 

a  =  {.80  +  .Ol(i-l)  i=l,...,20} 

K  =  220  =  1048576 

Lower  bound  based  on  3  edge-disjoint  paths 
Upper  bound  based  on  5  edge-disjoint  cutsets 


Fig.  2  Variance  Ratios  for  Alternative  Estimators 


•  1  v*r  iR(q)/v*r  fcbKU.p) 

▲  2  v«r  ^(qj/var  fcbKU.P) 

□  3  v«r  lK(q,p,9*Cq,p))/v*r  gbK(q,p) 


p  ■  .80 

K  -  1048576 


Appendix 


Proof  of  Lemma  1.  Observe  that 


E[#X)R(X,k,q,p)l  =  E  #x)R(x,k,q,p)Q(x,k,p) 

xt& 


Since  ^x)^(x)  =  ^(x)  and  ^x)^(x)  =  *(*).  one  has 


Also, 


E(^X)R(X,k,q,p)]  =  [g(qhgf(q)]/A(p). 


E[^(X)R(X^q,p)]2  =  E  ^(x)R2(x,k,q,p)Q(x,k,p) 

xejr 


xH,W 


Observe  that 


PJ(xtk,q)/P(x,k,p)  =  n  (^/Pj)*1  [( 1-qj  J’/O-Pi  )]ki  X\ 

i=l  1  1  1  1 


=  c(q,p){  n  (q?/cipi)Xi[(l-qi)2/ci(l-pi)]  1 


k.-x. 


80  that  {P2(x,k,q)/c(q,p)P(x,k,p),  x€JS}  is  a  p.m.f.  Expression  (12)  follows. 


Proof  of  T^nvma  2.  We  restrict  z  to  [gL(q),  g^(q)].  Consider  the  case  p>q  which  implies 
that  q*<q  so  that  g(q*)  <  g(q)  and  gL(q*)  <  gL(q).  In  this  case  h^z)  <  h2(z)  so  that  (40b) 

gives  the  tightest  upper  bound. 

Since  p<q  implies  q<q*  so  that  g(q)  <  g(q*)  and  g.(q)  <  g.(q*)  for  je{L,U},  one  has 
g(q  )  *  )>  «(q)]-  Also,  either  gL(q)  <  gL(q*)  <  g(q)  or  gL(q*)  >  gy(q).  Since  g  (q*) 


I  gj(<l)  implies  gjq  )  >  g(q),  h^z)  <  hj(?)  so  that  (40b)  gives  the  tightest  bound.  if  «l(q) 
<  %($*)  <  gu(q),  it  is  not  dear  whether  g(q)  <  gjq*)  or  g(q)  >  gL(q*).  Therefore,  (40a) 

gives  the  best  bound. 

Proof  of  Theorem  3.  The  function  hj  has  its  unrestricted  maximum  at  z  =  gy(q). 

Therefore, 


gt(^M<ohlW“hl(l°(q)) 


»f  gjq  )^[gL(q),gc(q)] 


ifsL(q  )e[gl(q),gu(q)]. 


The  function  hj  is  concave  with  its  maximum  at  *  <g^(q).  Therefore, 


.  35  ) 

gL(q  )<*<gp(q) 

Then  Theorem  •*  follows  directly  from  Lemma  2. 


if  gL(q  )<z  <gg(q) 


*  .  * 


if  z  <gL(q  )• 


Theorem  4.  Define  Wj(z)  =  \(z)/z2  and  w2(z)  =  h2(z)/z2  for  and  h2  as  in  (39) 
for  —«#<z<a#.  Let  ij  *  gg(q)  —  c(q,p)A(p)A(q*)/gg(q),  z2  =  [gJ(q)-c(q,p)A(p)A(q*)]/ 
[gfl(qH(q,P)]A(p)/2]  and  b2  =  2g(J(q)  -c(q,p)A(p).  Then 


^b(';V(q,P)  =  max  w  (z)  if  gL(q*)*[gL(q),g0(q)] 

gJqJSz^q)  h  u 


max  *  w  (z),  *  max  w„(z) 

gJqhzSg^q )  gL(q )  i  z<g^,(q) 


if  gL(q  )e[gL(q),gu(q)] 


i. 

,  ,  *i(z)  =  wi(8l(<0) 

*,,(<!)  S  •<*,(<!)  1  in 

=  *,(*,) 

if  zt<gL(q) 

if  z^gjq), 

u. 

max  *  w  (z)  =  w  (q)) 
%,(*)*«*%(«  )  1  1  L 

if  z^gL(q) 

■  "fa)) 

if  g^q^Zj&Jq*) 

and 

=  w1(gL(q*)) 

if  z^gjq*) 

Ui. 

,  •.“*  ,  "fa  = 

«i,(q  )<*<gD(q) 

if  z2<0  and  b2>0, 

*  w2(g0(q)) 

if  z2<0  and  b2<0 

=  w2(«i(q*)) 

if  0<Zj<gL(q*)  and  b2>0 

-  w2(*2) 

if  gL(q*)<z2<gu(q)  and  b2>0 

=  WjCgpCq)) 

if  z2>g^(q)  and  b2>0 

=  w2(gg(q)) 

* 

if  0<Zj<gL(q  )  and  b2<0 

=  max[w2(gL(q*)),w2(gD(q))] 

if  gL(q*)<z2<gu(q)  and  b2<0 
=  w2(gL(q  ))  if  ^  b2<0- 


Proof  of  ports  i  and  li.  Since 


—44— 

=  V*2  +  V  *-l 

where 

*  c(q,p)A(p)A(q*)  -  gjj(q)  and  bj  =  2gu(q), 

then 

*  — (2a1/*-Kb1)/*2  =  -  b1(-*1/z+l)/z2 

and 

dV  ,  , 

— 2(3a1/z+b1)/zs  =  2b.(-3z  /2z+l)/z3. 
dz4  *  1  11 

If  z1<0,  then  Wj,  is  convex  and  decreasing  on  [0,*)  and  w1  has  its  maximum  at 
a-^Cq).  If  tj>0,  then  w1  is  concave  on  [0,3z1/2]  and  z1<gQ(q)  so  that  w1  has  its 
maximum  on  (gjq),  gp(q)]  at  z  =  tl  if  z^g^q),  g^q)]  and  at  z  =  gw(q)  otherwise. 
The  result  for  part  ii  follows  immediately. 

Proof  of  in.  Since 

w2(«)  *  hj(«)/z2  =  a2/z2  +  b2/z  - 1 

where 

*2  *  c(q,p)A(p)  g0(q*)  -  gj(q)  and  b2  =  2g0(q)  -  c(q,p)A(p), 

then 

dw,  0 

~t)2(-z2/z+l)/z 

and 

d2w 

2b2(-3z2/z+l)/z3. 

Consider  the  interval  [g^q*),  g^(q)].  If  z2<0  and  b2>0  then  w2  is  convex  on  [0,®)  and  its 
maximum  occurs  at  z  =  g^(q  ).  If  z2<0  and  b2<0  then  w2  is  concave  on  [0,®)  and  has  its 
maximum  at  z  *  gg(q).  If  z2>0  and  b2>0,  then  w2  is  concave  on  [0,3z2]  so  that  the 


maximum  occurs  at  *  »  gjq  )  if  SjSg^q  ),  *  *  *2  if  gL(q  )<^<gu(q)  and  at  z  =  g^q)  if 
Sj^g^q).  If  *j>0  and  b2<0,  the  maximum  occurs  at  z  =  gg(q)  if  z2<gL(q*)>  at  2  =  gL(q*)  if 
a^(q)  and  at  a  =  maxfw^g^q*)),  w2(g0(q))]  if  gL(q*)<z2<gy(q). 


Theorem  5.  Let  w3(z)  *  h1(*)/(l-z)  and  w4(z)  =  h2/(l-z)  for  hL  and  h2  as  defined 

i&  (39).  Let 


*j  *  g^Q)  +  c(q,p)A(p)  A(q  )/[l— g^Cq)), 


*.  *  1 


Then 


-2{c(q,p)A(p)(i-gD(q  )]  +  (l-gD(q)]2}/{c(q,p)A(p)  +  2[l-gu(q)]}. 


^b(<W»)  <  w**(q,p)  =  max  w  (z)  if  gL(q*)*[gL(q),g0(q)] 

gL(q)^*ig0(q) 


where 


=  max  max  * « 

[gL(q)<z<gL(q  ) 

r  (z)  *  max  w  (z) 

s  )<'Sg„(<i) 

■f  sl(q*)E|gl(q)^u(q)] 

max  ws(z)  =  w  (g  (q)) 

=  w^q)) 

ifz3<i 

,  xmAX  ,  *VWS(*)  =  w*(8l^*)) 
gl(q)<«<gL(q)  3 

if  Z3~1 
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*  max  w  (s)  =  w  (g  (q  )) 

^(1  )<«<%,(<!)  4  41 


if  *4<gL(<l ) 


“  w4<*4) 


if  gL(q  )<z4<gu(q) 
if  z4>g0(q)- 


whore 

and 

then 

and 


5.  Since 

wj(»)  =  +  bj/(l-z)  ~  1 

03  *  c(q,p)A(p)  A(q*)  -  (l-g^q)]2 

bj  =  2(1  g|j(q)], 

d--^  [w*H 


dz  (l-*)41 


ST-"^  [3Vb»(1-*)+1]  =S?h(I^»)/2(1^)+1]' 


Observe  that  *s  =  1  +  2aj/b3  >  gjj(q).  if  *3>1,  then  w3  is  convex  on  (-»,1],  having  its 
maximum  on  (gjq),  gu(q)]  and  on  (gL(q),  gL(q*)]  at  2  =  gL(q).  If  z3<l,  then  w3  is  convex 
on  ((3*3—1  )/2,l]  and  wx  has  its  maximum  on  (gjq),  g^q)]  at  z  =  g^q)  and  on  (gL(q), 
gL(q  )]  at  2  =  gjq  ),  establishing  i  and  ii. 

Proof  of  Ui.  Since 

w4(»)  =  +  V^-*)  ~ 1 

where 

*4  =  cfqrfOAfrHl-g^q*)]  +  (l-gu(q)]2 

and 


b4  =  2(1— g^q)]  +  c(q,p)A(p), 


m  ..  t 


[2a4/b4(l-»)+lJ  =  [-(l-*4)/(l-«)+lJ 


dV  2b 


d?2  =  (l^j3  [3Vbl*1_z)+2]  =  [-3(l-*4)/2(l-z)+l]  • 


(1-*)“ 


Consider  the  interval  [gL(q  ),  gg(q)].  Since  z4<l  w2  is  concave  on  [ — (1— z4) /2,1].  Then 
w4  has  its  maximum  at  z  =  gL(q*)  if  z4<gL(q*),  at  z  =  z4  if  gL(q*)<z2<gu(q)  and  at  z  = 

%(<!)  ^  *4>g0(q)- 

Proof  of  Theorem  8.  Observe  that 


R*(x£,q,p)  =  n  (qi/Pj)  1  ((l-q.)/(l-p.)] 
i6  Jr 


k.-x. 
i  i 


has  the  alternative  form 


R*(x,k,q,p)  =  [  n  (q./p.)  l]  exp[-  E  a  £  (1— y  )] 
i€<*  1  1  i6<^  1  leaf.  J 


E  y.=x 


The  condition  ^(x)  =  0  requires  that 

^ 2 1 

and  the  condition  $x)  =  1  requires  that 
E  (1-yJ  <  |tf|  - 1 


ieJP*. 


V 


V  #  e jr. 
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Since 


max  exp[—  E  a.  E  (1— y .)]  =  exp[—  min 
y  i  e^T*  j€*r  1  y 


E  a. 
i  ejr* 1 


j€*. 


d-y,)], 


one  has  the  integer  program  (47)  with  y^  =  (l-zp. 


